First, we load the data into the workspace.

Following is a summary of the dataset.

1. Exploratory analysis

Temporal: Frequency domain

The analysis in the frequency domain reveals a strong harmonics with a period of 24 hours. Therefore, it is reasonable to

freq.plot(checkin.NY,title="New York City")

plot of chunk unnamed-chunk-4

Temproal: Time domain

The distribution of the check-in categories accross one period (24 hours).

time.distribution.plot(checkin.NY,title="New York City")

plot of chunk unnamed-chunk-5

Temporal: Radial plot

The distribution of the check-in categories in 24 hours in radial plots.

time.radial.plot(checkin.NY)

plot of chunk unnamed-chunk-6

Spatial

have a look at a intuitive spatiotemporal visualizaiton of the data…

basemap = map.plot("../../global/data/shapefiles", 
                   "NYC_borough_boundaries_WGS84",
                   alpha=0.1,size=0.3,color="grey")

saveGIF({
  point.animation.plot(checkin.NY,title="New York City", 
                       basemap=basemap,axis.size=18,
                       title.size=24,legend.size=18,plot.size=3)
}, interval = 0.5, movie.name = "spatiotemporal.gif", 
ani.width = 1500, ani.height = 1200)

What will the statistics say when we neglect the temporal factors? First of all, if we try to find spatial clusters based on differnt checkin categories, the distribution of the founded clusters are quite differnt. It indicates that each category has differnt correlations with the geographic space.

if(!exists("basemap")){
    basemap <- map.plot("../../global/data/shapefiles", 
                   "NYC_borough_boundaries_WGS84",
                   alpha=0.1,size=0.3,color="grey")
}
## OGR data source with driver: ESRI Shapefile 
## Source: "../../global/data/shapefiles", layer: "NYC_borough_boundaries_WGS84"
## with 5 features and 4 fields
## Feature type: wkbMultiPolygon with 2 dimensions
plots <- lapply(split(checkin.NY,checkin.NY$cate_l1),function(ci){
    # find out clusters for the type of category
    clusters = spatial.clustering(ci)
    
    centers = clusters[["centers"]]
    points = clusters[["point.unique"]]
    wss = clusters[["wss"]]
    pct = clusters[["pct"]]
    lbl = paste(nrow(centers),"clusters;\nWSS:",
                  formatC(wss,digits=2,format="f"),
                  "(",format.percent(pct),")")
    
    # add basic points
    gg.map <- point.plot(points,x="lon.x",y="lat.x",alpha=0.05, basemap=basemap)
    
    # add cluster information
    gg.map <- point.plot(centers, x="lon.center",y="lat.center",
                         color= "cid.ordered", color.concrete=FALSE,
                         point.size = log(centers$size,5),alpha = 0.7,
                         basemap=gg.map,
                         xlim=range(checkin.NY$lon,finite=TRUE),
                         ylim=range(checkin.NY$lat,finite=TRUE))
    
    # some plot configuration
    gg.map <- gg.map + ggtitle(ci[1,"cate_l1"]) +
        theme(legend.position="none") + 
        geom_text(aes(x = -74.13, y = 40.87), label = lbl, size=2) 
    
})
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
# png("categorized_clustering.png",width=15*ppi,height=6*ppi,res=ppi)
grid.arrange(plots[[1]],plots[[2]],plots[[3]],plots[[4]],
             plots[[5]],plots[[6]],plots[[7]],plots[[8]], 
             plots[[9]],plots[[10]],nrow=2, ncol=5)

plot of chunk unnamed-chunk-8

# dev.off()

We could also see the most domimant categories in each part of the entire city.

# prepare polygon data
SPDF = readOGR(dsn = "../../global/data/shapefiles", layer = "NYC_zipcode")
## OGR data source with driver: ESRI Shapefile 
## Source: "../../global/data/shapefiles", layer: "NYC_zipcode"
## with 236 features and 6 fields
## Feature type: wkbPolygon with 2 dimensions
# intersect the checkin data with the each postal code region
checkin.NY = point.in.poly(checkin.NY, SPDF, copy.attr="POSTAL")
## [1] "Spatial data has been prepared."
## [1] "The information from the polygon has been assigned to the points."
# get the distribution of categories in each postal code region
cate.distr=cate.distr.in.poly(checkin.NY)
## [1] "The distribution of category has been assigned to each polygon."
# and use that information to implement the SPDF
SPDF = poly.with.category(SPDF, cate.distr=cate.distr, poly.attr="POSTAL")

# plot

mapdf=df.from.spdf(SPDF)
# mapdf$density=apply(mapdf,1,function(i){i[i["cate.dom"]]})
# mapdf$density=as.numeric(formatC(mapdf$density,digits=1,format = "f"))
png("main_category.png",height=6*ppi, width=16*ppi,res=ppi)
grid.arrange(
map.plot(mapdf=mapdf,
         more.aes=aes_string(fill="Category.1st"),
         color="grey",size=0.3,alpha=0.9)+
    theme_bw()+
    xlab("")+ylab("")+
    theme(legend.title = element_text(size=8),
          legend.text = element_text(size = 8)),
map.plot(mapdf=mapdf,
         more.aes=aes_string(fill="Category.2nd"),
         color="grey",size=0.3,alpha=0.7)+
    theme_bw()+
    xlab("")+ylab("")+
    theme(legend.title = element_text(size=8),
          legend.text = element_text(size = 8)),
nrow=1, ncol=2, widths=c(1,1) )
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
dev.off()
## pdf 
##   2

Meteorological: Correspondence analysis

cate.conds = xtabs(~conds+cate_l1, data=checkin.NY)
#prop.table(cate.conds, 1) # row percentages
#prop.table(cate.conds, 2) # column percentages
fit <- ca(cate.conds)
#print(fit) # basic results
summary(fit) # extended results
## 
## Principal inertias (eigenvalues):
## 
##  dim    value      %   cum%   scree plot               
##  1      0.003244  52.2  52.2  *************************
##  2      0.002082  33.5  85.8  ****************         
##  3      0.000511   8.2  94.0  ****                     
##  4      0.000147   2.4  96.4  *                        
##  5      7.6e-050   1.2  97.6  *                        
##  6      6.7e-050   1.1  98.7                           
##  7      5.1e-050   0.8  99.5                           
##  8      2.4e-050   0.4  99.9                           
##  9      8e-06000   0.1 100.0                           
##         -------- -----                                 
##  Total: 0.006210 100.0                                 
## 
## 
## Rows:
##      name   mass  qlt  inr    k=1 cor ctr    k=2 cor ctr  
## 1  | Cler |  501  779   88 |   22 452  76 |  -19 327  86 |
## 2  |  Fog |    2  634   31 | -265 628  37 |  -26   6   1 |
## 3  | Haze |   22  975  178 | -200 798 272 |   94 177  94 |
## 4  | HvyR |    5  917   75 |  208 421  61 |  226 496 111 |
## 5  | HvyS |    0  116    5 | -170  20   0 |  377  96   1 |
## 6  | LgFR |    0  714   29 | -590 712  40 |   28   2   0 |
## 7  | LghR |   33  947   88 |   77 362  61 |   98 586 154 |
## 8  | LghS |    7  246   21 |  -62 203   8 |   29  44   3 |
## 9  | MstC |   80  738   38 |  -31 319  23 |  -35 420  47 |
## 10 | Ovrc |  232  957  107 |  -25 217  45 |   46 740 237 |
## 11 | PrtC |   55  592   41 |   38 315  25 |  -36 277  34 |
## 12 | Rain |   10  943   67 |  163 632  81 |  114 312  62 |
## 13 | SctC |   49  874   97 |  -59 287  53 |  -85 587 169 |
## 14 | Snow |    4  850  134 | -421 849 218 |    4   0   0 |
## 
## Columns:
##      name   mass  qlt  inr    k=1 cor ctr    k=2 cor ctr  
## 1  | ArtE |  105  763   76 |   49 526  76 |   33 237  54 |
## 2  | CllU |   17  677   27 |  -75 580  30 |   31  97   8 |
## 3  | Evnt |    5  557   20 |  -24  23   1 | -113 534  32 |
## 4  | Food |  294  688   48 |   18 301  28 |  -20 388  56 |
## 5  | NghS |  122  955  206 |   93 826 325 |   37 129  80 |
## 6  | OtdR |   69  890  200 |   -9   4   2 | -127 886 528 |
## 7  | PrOP |   98  986  272 | -123 878 457 |   43 108  87 |
## 8  | Rsdn |   27  352   42 |   23  54   4 |   54 298  38 |
## 9  | ShpS |  156  661   49 |  -29 430  41 |  -21 230  34 |
## 10 | TrvT |  107  791   59 |  -33 318  36 |   40 473  84 |
#plot(fit) # symmetric map
plot(fit, mass = TRUE, contrib = "absolute", map ="rowgreen", 
     arrows = c(TRUE, FALSE)) # asymmetric map

plot of chunk unnamed-chunk-10

3. Model - derivation and corresponding functions

Under the assumption \(H=i\) is independent from \(W=j\) \[P(C=k|H=i,W=j)=\frac{P(C=k,H=i,W=j)}{P(H=i,W=j)}=\frac{P(H=i,W=j|C=k)*P(C=k)}{P(H=i)*P(W=j)} (1) \]

since \(H=i\) is independent from \(W=j\), \[Exp[P(H=i,W=j|C=k)]=P(H=i|C=k)*P(W=j|C=k) (2)\]

therefore, \[Exp[P(C=k|H=i,W=j)]=Exp[ \frac{P(H=i,W=j|C=k)*P(C=k)} {P(H=i)*P(W=j)}] \\\ =\frac{P(H=i|C=k)*P(W=j|C=k)*P(C=k)}{P(H=i)*P(W=j)} \\\ =\frac{\frac{P(H=i,C=k)}{P(C=k)}*\frac{P(W=j,C=k)}{P(C=k)}*P(C=k)}{P(H=i)*P(W=j)} \\\ =\frac{P(C=k|H=i)*P(H=i)*P(C=k|W=j)*P(W=j)}{P(H=i)*P(W=j)*P(C=k)} \\\ =\frac{P(C=k|H=i)*P(C=k|W=j)}{P(C=k)} (3)\]

\[P_{u}(C=k|H=i)=\frac{\Phi_{u}(C=k,H=i)}{\Phi_{u} (H=i)} (4)\]

get.temporal.impact <- function(dataframe,hour){
    # dataframe.in.hour = checkin.single[which(checkin.single$hour==hour),]
    dataframe.in.hour = dataframe[which(dataframe$hour==hour),]
    phi.h = nrow(dataframe.in.hour)
    
    list.category = split(dataframe.in.hour, dataframe.in.hour$cate_l2)
    sapply(list.category, function(i){
        nrow(i)/phi.h
    })
}

\[P_{u}(C=k|W=j)=\frac{Intercept(C=k,W=j)}{\sum Intercept(C,W=j)} (5)\]

get.meteorological.impact <- function(fit,conds){
    
    conds.id = which(fit[["rownames"]]==conds)
    ref.vec = fit[["rowcoord"]][conds.id,1:8]
    cate.all = fit[["colcoord"]][,1:8]
    
    intercepts = apply(cate.all, 1, function(x){ 
        (x[1]*ref.vec[1] + x[2]*ref.vec[2] + x[3]*ref.vec[3] + x[4]*ref.vec[4] +
             x[5]*ref.vec[5] + x[6]*ref.vec[6] + x[7]*ref.vec[7] + x[8]*ref.vec[8] ) /
            (ref.vec[1]^2 + ref.vec[2]^2 + ref.vec[3]^2 + ref.vec[4]^2 +
                 ref.vec[5]^2 + ref.vec[6]^2 + ref.vec[7]^2 + ref.vec[8]^2 ) 
        } )
    
#     vec = intercepts / sum(intercepts)  !!!wrong!!! intercetp can be negative
    vec = (intercepts - min(intercepts)) / (max(intercepts)-min(intercepts)) # scale into [0,1]
    names(vec) = fit[["colnames"]]
    
    vec
}

get.meteorological.impact2 <- function(dataframe,conds){
    
    dataframe.in.conds = dataframe[which(dataframe$conds==conds),]
    phi.c = nrow(dataframe.in.conds)
    
    list.category = split(dataframe.in.conds, dataframe.in.conds$cate_l2)
    sapply(list.category, function(i){
        nrow(i)/phi.c
    })
}

\[P_{u}^{*} (C=k|W=j)= w_{j}*[P_{u}(C=k|W=j)-\bar P_{u}]+\bar P_{u}\]

get.weather.weight <- function(fit){
    
    conds.all = fit[["rowcoord"]][,1:2]
    
    mag = apply(conds.all, 1, function(x){ 
        sqrt( (x[1]^2+x[2]^2) )
        } )
    
    mag / sum(mag)
}

get.weighted.meteorological.impact <- function(fit, conds, weight){
    
#     cate.conds = xtabs(~conds+cate_l2, data=dataframe)
#     fit <- ca(cate.conds)
    
    unweighted = get.meteorological.impact(fit,conds)
#     weights = get.weather.weight(fit)
    
#     conds.id = which(fit[["rownames"]]==conds)
    
#     vec = weights[conds.id] * (unweighted-mean(unweighted)) + mean(unweighted)
    vec = weight * (unweighted-mean(unweighted)) + mean(unweighted)
    names(vec) = fit[["colnames"]]
    
    vec
    
}


get.weighted.meteorological.impact2 <- function(dataframe, conds, weight){
    
    
    unweighted = get.meteorological.impact2(dataframe, conds)

    weight * (unweighted-mean(unweighted)) + mean(unweighted)
    
    
}

\[P(C=k)=\frac{\Phi_{u} (C=k) }{\Phi_{u} }\]

get.denominator <- function(dataframe){
    
    phi.h = nrow(dataframe)
    
    list.category = split(dataframe, dataframe$cate_l2)
    denominator = sapply(list.category, function(i){
        nrow(i)/phi.h
    })
    
    denominator

    
}

\[E[P(C=k|H=i,W=j)]=\frac{P(C=k|H=i)*P(C=k|W=j)}{P(C=k)}\]

update @2014.09.04

update @2014.09.05